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SUMMARY 


An algorithm and ILLIAC computer program, developed for the simulation 
of homogeneous incompressible turbulence in the presence of an applied mean 
strain, are described. The turbulence field is represented spatially by 
a truncated triple Fourier series (spectral method) and followed in time 
using a fourth-order Runge-Kutta algorithm. Several transformations are 
applied to the numerical problem to enhance the basic algorithm. These 
include 

1. Transformation of variables suggested by Taylor’s sudden-distortion 
theory 

2. Implicit viscous diffusion by use of an integrating factor 

3. Implicit pressure calculation suggested by Taylor’s sudden- 
distortion theory 

4. Inexpensive control of aliasing by random and phased coordinate 
shifts 


INTRODUCTION 


The primary difficulty in the numerical simulation of homogeneous 
turbulence is that the nonlinearity of the equations of fluid motion excites 
a large range of scales (i.e., a large ratio of largest to smallest scale) 
of motion in both space and time. The computer resource required for a 
complete simulation is proportional to the product, over all space-time 
dimensions, of the range of computed scales of each dimension. These scale 
ranges increase with Reynolds number (R) , and their product increases so 
rapidly, in three space dimensions, that only the weakest experimentally 
studied turbulence can be simulated completely on today’s computers. 

The overall range of scales continues to increase indefinitely with 
Reynolds number. However, at a sufficiently high Reynolds number the scales 
of motion can be grouped, in order of decreasing scale, into three distinct 
ranges: the energy-containing range, the ’’inertial” range, and the dissipa- 

tion range (fig* 1). Further increases in Reynolds number Increase only the 
inertial range. The range of energy-containing scales, which determine the 
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Ic'aturcs ot turbulence of engineering interest, is bounded «s R ®. and the 
motion in these scales becomes independent of tl.e motion at smaller scales 
At somewhat lower U, the inertial and dissipation ranges merge, but still'do 
net alfect the energy-containing range. At sufficiently low R, dlsslpatJoj 
occurs In the energy-cemtaining range Itself. This physical description 
of the scale dependence upon Reynolds number is encouraging because^lt Indi- 

L inci^ded " energy-containing scales of motion nLd to 

is ^ ^ high Reynolds number turbulence simulation. The difficulty 

terms lA ‘ all the scales are coupled through the nonlinear 

terms in the governing equations and, although we know that physically 

iman'er^c"* energy-containing range is uncoupled from the 

smaller scales, we do not know how to uncouple it mathematically. 

an-fo J^ange of statistically interdependent scales increases with the 

anisotrooL 1^ i engineering Interest are 

^ important to determine the nature and magnitude of the 
additional computational difficulty posed by anisotropy. 


THE NUMERICAL SIMULATION 


NavleJ-StoZ‘‘cL'i°;f praaantad here is an unsteady Incompressible 

Navier Stokes code that runs on the ILLIAC IV computer. The program co-nmii-,.., 

e evo ut on in time from an arbitrary homogeneous turbulence field in 
the presence of a simple class of spatially-linear mean flows! tL slnula- 
Hff ^ spectral decomposition similar to that of Orszag (ref. 1) but 
differing in detail. The primary purpose of this report is to present the 
simulation algorithm in detail sufficient to allow its use by others III 
lencrf presented to study weak (low Reynolds^number) 'turbu- 

coZr.M """ presented. The magnitude of the 

P on ( 3.g, 2) requires a computer at least as fast as a CDC-7600. 


THE EQUATIONS OF MOTION 


The equations governing the flow of a viscous 
the familiar Navier-Stokes equations 
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where (u,v,w) Is the vclnclty vector , p Is the pressurt- -deiisi t v ratio, 

V is the kinematic viscoalty, and Hubsi.ri))ts denote di iferent lat i on . 

We wish to simulate numerically the el feet of a simple class of ii,i,.()‘;ed 
strains on a homogeneous field of turbulence. The strain field is given by 


(u.v.w) = [xa(t),yb(t),zc(t)l 

where a + b + c = 0 as required by continuity. It is convenient to intro- 
duce the following transformation of the dependent variables: 


u “ ax + A^'^^u 
V = by + B^^^v 
w = cz + C^^^w 


^ ( 2 ) 


-i[( 


da 

dt 


+ a 




Ik 

dt 


+ b 


* ( 


dc 

dt 


+ c 


•y] 


+ P 


where a(t), b(t), and c(t) are the arbitrary time-dependent strain rates 
Imposed, and the resulting Inverse square strains are 



a dt 


b dt 


> (3) 


~2 f ^ c dt 

C(t) = e 

y 

It follows from the continuity condition that material volumes are invariant, 
(i.e., ABC = 1). Explicit spatial dependence of the resulting system of 
equations is eliminated by the following transformation of independent 
variables : 


X = A^^^x 
y as B^^^y 

2 » C^^^z 


( 4 ) 


The equations of motion for the transformed turbulence field arc then 



Uj. + A(uu)‘ + B(vu)^ + C(wu).; + P' = 


V 4- A(uv)- f B(vv)~ + C(wv); + p- = v(Av-- + Bv— + Cv--) 
t X y 2 ' y XX yy zz 


Wj. + A(uw)^ + B(w)' C(wv).^ + p; « V(AW' ' -i- BW • + CW'") 

t X y z ' 2 XX yv zz 


Au;' + Bv' + Cw- = 0 
X y z 
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Tc is interesting to note that Taylor’s theory (see ref. 2) for the 
sudden distortion of a turbulent field is particularly simple when expressed 
in these transformed variables, because strain mte does not appear. Suppose 
that a field of turbulence (state 1) is rapidly strained. The resulting 
field (state 2) is defined as 


fi<2) . o«) ^ r c,^ 


dt , etc. 


where t is a time characteristic of the imposed straining period. 
Continuity must be satisfied throughout the straining process and, in 
particular, in the initial and final states of strain. Hence, if the total 
strain imposed is given by its inverses A,B,C, we have 

i<‘> + JP) + i!') .0 at t 

X y z 

Au|^^ + =•- 0 at t + t 

X y Z 

As the strain becomes sudden (t 0) A,B,C, u, v, and w remain bounded but 
their time derivatives do not. During the straining process the momentum 
equations (5) degenerate, as t U, to 

Ot + p. . 0 

+ pp . 0 

“t + Pj * 0 

SO that the velocity jump in the transformed variables is irrotatlonal , 
that is 
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and similarly 
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Thh contirmity condition then determines the potentlnJ as the solution 

=■ + BV'*^ + 

XX yy zz x y 
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which completes the solution for the veiority field at state 2, 


The above transformations seem to be the natural ones for the study of 
the effect of uniform imposed strain on a homogeneous turbulent field > 
regardless of the rate at which the strain is imposed. A more general set of 
transformations can be used when the mean strain-rate matrix is not diagonal, 
and also when the mean vortlclty is nonzero. 


NUMERICAL APPROXIMATIONS 


We wish to simulate a spatially homogeneous turbulence field in an 
infinite space, and this suggests that we represent the field spatially as a 
Fourier series. The resulting field is periodic in all three space 
dimensions, with correspondingly periodic spatial correlations. However, 
if these correlations decay to negligible magnitude within the period, 

(c.g., if the integral scale is much smaller than half the period), the error 
du3 to the finite period should be small. In practice this requirement is 
difficult to satisfy with the resolution allowed by today’s computers. 


In this section we develop the equations in more detail and describe 
the integration process as programmed. Let = uu, ^13 * 

etc., tilde (-) denote the three-dimensional Fourier transform, and 

^2* ^3 wave numbers in the x, y, and z directions, respectively. 
The equations (5) in wave space are then 


“t 
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IkjAXjj + 


ikjCtj, H- 

Ikjp 



ikjAtj2 + 

ik2B^22 + 

lk3ri^3 + 

ik2P 


-v(Akj2 + 

-v(Ak/ + Bk^^ 


+ Ckg^)u 
+ Ck3^)v 


V(6) 


Wj. + ikjAtjj + “ -vCAk^^ + Bk,^ + Ckj‘^)w 

ik Au + Ik^Bv + ik Cw = 0 

1-3 J 

The linear terms are combined by multiplying the equations by the integrating 
factor 
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^ H dl vk.,^ C 

F(k.t) - f ‘' e ' '’ c 


giving 

(Fu) + I'fikjATjj i- ik^Btj^ + Ik^Ctj.^ + ik^p} = o 

~ (Fv) + F{lk,A^j 2 + ^^2^^ 22 3 + 

d 

^ (FV) + F{ikjATj3 + ik2BT~23 + ik3C^3^ + ik3p) = 0 
ikjA(Fu) + ik 2 B(Fv) + ik 3 C(Fw) = 0 

Now multiply the first equation by ik^ , the second by ik etc., to obtain 
(let U « ikju, V = lk 2 V, W = lk 3 w) ^ 


~ (FV) - F{kjk2A^^^ + k^2B,-^^ + k2k3C^23} + k^^Fp 

(9; 

dt + k 32 pp 

A(FU) + B(FV) + C(FW) = 0 

(The purpose of this transformation of dependent variables is discussed 
later on; note that ^i»^2’^3 ~ ^ special cases.) 

The usual procedure for the computation of p requires the time 
differentiation of the continuity condition. However, we want the algorithm 
to handle impulsive strains correctly (Jumps In A, B, and C), that is 
according to Taylor's sudden distortion tfieory, so we need to avoid the 
differentiation. We thus define a potential ^ as 
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Fp dt 


and absorb It into the time-advanced varld]>les. 
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and the continuity condition becomes 

; = F-l/ Ai+_BYJ_CZ__ \ 

Uk ^2 + bi ^^2 4. J 

where 

F~^X = U + kj2^ 

F“1y = V + 

F'^Z = W + k^2j, 

The t's are functions of u, v, w only, so that, if u, v, w are known 
at the beginning of a time step and satisfy the continuity condition, we may 
advance X, Y, Z, However, to form time derivatives of the advanced X, Y, Z 
we must form from (12) advanced values of u, v, w, and this requires the solu- 
tion (11) for (P at the advanced time. This is done using the continuity 
condition at the advanced time, and does not require its time differentiation. 

At the beginning of a time step t = 0, and we have 

F=1,^ = 0,X=U,Y = V,Z = W 

The equations for X, Y, and Z are integrated over the time step, and the 
final values are used in equations (11) and (12) to produce final values of 
U, V, and W. The origin of time; is then shifted to the final time giving the 
proper initialization for the next time step. 

Spatial dif fe^ntiatlon is a point operator in wave space but multiplica- 
tion (e.g., Tj 2 * uv) is not, and the most efficient means of forming the 
Fourier transform of a product from the transforms of its terms is to return 
to physical space by inverting the transforms, form the product, <»nd then 
transform the result back to wave space. Unfortunately, the transformation 
of the product back to wave space Introduces an error due to spectral 
truncation. 
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Tlio truncation crrorn arc* ;r,oHl oaHlly cliMnonsl rat cd in one .spaLiai 
dlmonsLon. The* ropmsmt/it I f>n «>r ilu* product of two Fourier .sorien a, h 
(In c^oniplox rorm) as a Fourlt*r seric's c U\ p/iven hy tlic* (inrinitc) 
(U)nvoliitlon sum 






a, b 
k-s s 


However, the process of InvertiiiR finite transforms a and b, forming the 
product ab, and then taking its finite transform results instead In two 
sums: 



, b 
k-s s 


^ ^ktM-s^’s 
s 


The first sum represents a contribution (incomplete due to truncation) to ab 
correctly attributed to wave number k. The second sum also represents a 
contribution to ab, but it is actually a contribution not to k, but to 
k±M, wave numbers beyond those allowed by the length (M) of the finite 
transforms used. This is the "aliasing" error. Now it may be argued that 
because aliasing errors do not account for all of the truncation error, 
suppression of the aliasing error is not cost effective so far as accuracy is 
concerned. However, in the algorithm used here, the aliased terms can lead 
to nonlinear instability, and their control is essential. 


Now consider the effect of a shift of the physical coordinate system. 

In wave space this amounts to multiplication by where -A is the 

amount of coordinate shift. If we use e"^*^ to shift aj^, bp prior to 
inverting them to physical space, form the product ab on the shifted grid, 
transform back to wave space, and finally shift coordinates back with 
we obtain 


■E 


a, b + e 
k-s s 


tlMA 


' ^k±M-s^s 


The first (alias-free) sum is invariant under these shifts, but the -jecond 
sum, the aliased one which we wish to suppress, has a phase dependency on A 
and can be eliminated. For example, if two evaluations are made, one with 

o * 1 and the other with e“^^ = -1, the alias-free result is one-half 
their sum. The second sum (which is multiplied by the phase factor) itself 
vanishes identically for |k| < N, (N < M/1) if modes of a and b outside of 
this range are nulled prior to inversion, and transforms of length M aie 
retained. Thus two independent procedures are available for alias suppression. 

The extension of these procedures to three dimensions gives for each 
eight terms, seven of which represent aliasing errors. The aliased terms 
are tlassifled according to the number of dimensions in which aUasiii); has 

+ i knAn 

occurred. We have |k (k, , k ,, k^), 1 





(. 1 1 las-f j (H-) 

(k in}’ 1 y-nl 1 as<’d) 
(<l<ml)] y-;tl iriHod ) 
(t riply-.i l last’d ) 

All of the aliased sums (Sj , . . . S^) vanish If modes havln)- any > Ni 
are nulled. The doubly and triply aliased sums (5,^, • . S_) vanish if 

modes havinjs any two > Nj^ are nulled. The triply aliased sum (Sy) 
vanishes if modes having all thvoci k^ > are nulled. Alternatively one 
can evaluate the convolution eight times using the eight combinations of 
®x» ®y» ®z “ eliminate the aliased terms. Note that suppres- 

sion by the latter means requires eight evaluations to eliminate all of the 
aliased terms. One can also, as suggested by Orszag (ref. 1), remove 
S^, . . . Sy by truncation and the remaining single aliases by coordinate 
snift with two evaluations. We are faced with the choice between losing 
information (truncation) or losing computational speed (multiple evaluations). 

We have, following Orszag, eliminated doubly and triply aliased sums 
by truncation, though the truncation used here differs slightly from that of 
Orszag who nulls modes having ii * h. ^ 2(M/3)^. We have not exactly 
eliminated the remaining single aliases, due to the computational cost of 
the double evaluations required. Instead we have used the fact that the 
Runge-Kutta algorithm requires pairs of evaluations at each half step and 
that by using a shifted grid for the second evaluation we reduce the total 
alias error for the pair by a factor of At^. The possibility of nonlinear 
instability is further reduced by insuring that the Oj for the first 
evaluation in a pair are not correlated with those of other pairs. This is 
easily accomplished by the use of a unlform-random-number generator during 
computation of the phase factors. 


+ o.s, + n,H, t 0 

f 0,0y0.^Sy 


DATA MANAGEMENT 


In large simulations the high-speed random-access memory of the computer 
cannot hold the entire data base of the problem (in the present code it holds 
6% of it). In this case the high-speed memory may only be able to hold a few 
lines of the mesh (e.g., all values of kj for a few k 2 , kj values), and 
it is convenient to transform and take derivatives only along those 1 ines. 

In general, separate passes over the data base are required for each spatial 
dimension. The directional order in which operations are performed then 
determines the required number of passes over the data base. We will 
demonstrate how this number may be reduced in a spectral algorlilim. 

Consider the evaluation in wave space of (uv)x and (uv)", wliich is 
required in equation (5). The transforms of u and v are Inverted In the 
X, y, and z directions, each direction requiring a separate pass over the 
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data bast'. On t la* last (:/,) pass of t his .st*(jUi*nrf* wi* a lso l orm. In jjliysii'al 
s[>ac(', Lho uv )>rudu<‘t .md I lu-n t iansrurui back to wave spa« t' in r ho 
/, diroction. In prinrJpIo there remain only the x and y transforms and 
the* mult Ipllc/it ions by ik^ and 1ky to form the derivatives in tiie x and 
y directions. The problem la that, under our <’onst raints , transforms and 
derivatives can only be taken in the direction ef the* p,rid lines held In 
fast memory. Under these constraints we must cither perform three transforms 
and two di*r ivat J vos In two passes, or two transforms and two derivat Ives In 
three passes, 11 the constraint on the derivative la absent the results can 
be obtained in two transforms and two derivatives in two pasai*s. This 
constraint can be removed only if four lines of the mesh can be held simul- 
taneously in fast memory (so that all eiglit real numbers representing wave 
number k are present). The TLLIAC fast memory Is sufficiently large to 
accommodate four mesh lines, but not within a single processing element (PE), 
so that differentiation would require communication across the PE’s. We 
have instead used a slightly altered set of dependent variables that avoids 
this problem altogether. 

If the X momentum equation is differentiated with respect to x, and 
the y momentum equation with respect to y, the uv stress term appears as 
(uv)^y in both equations, and its evaluation under the constraints requires 
two transforms and two derivatives in two passes. Hut two extra integrations 
(of and v^) are then required to form u and v in physical space; 
however, since integration and differentiation cost far less than cither a 
transform or an I/O pass, this method is quite efficient. To avoid loss, 
upon differentiation, of information in a Fourier mode having a null wave 
number we simply do not multiply that mode by its wave number (l.e., zero) 
and similarly when we integrate it we do not divide by its wave number. What 
this amounts to is that, Instead of the usual spectral dependent variables 

u(k^ ,k^,k,) 
v(kj ,k.,,k J 
w(kj ,k^,k^) 

we use 


u (0 , kjj » k 2 ) , 

ik^u(kj ,k2,k^) , kj 

vCkj.O.kg) , 

ik^ V (k 1 , k^ , k ^ ) , 

w(kj,k^,0) , 

Ik ^w (k 1 , k2 , k *;^ ) , 


Use of these variables simplifies the continuity condition and minimizes the 
number of transforms and passes over the data base. 
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INITIAL CONr)TTION.S 


The Initial velocities are chosen randomly, subject to tlie constraints 
of continuity and a specified energy spectrum, in detail, tlie real and 
imaginary parts of the Fourier velocity amplitudes are selected ran- 

domly from a uniform distribution over the circle that is the Intersection 
of the sphere (having surface area proportional to 1’. (k)) determined by 
the desired energy spectrum and the plane (normal to k) determined by 
continuity. For example, consider the (real) spectral mode (kj,kj,,k 3 i 0) 


Uj^ = fj cos kjX sin kg',/ cos k^?. 
v “ cos k,x sin k,y cos k.z 

y \ j 

= fg cos kjX sin k^y cos k^x 

The algorithm described previously advances the vector _f in time given the 
initial values (k^.k^jk^ = 0 are special cases) 

. / kjk 

= c(kj2 + ~*1T 

/k k 

f^ = c(kj2 + k^^)“^/ 2 k 2 (-— sin i(; - k 
1/2 ^3 

= -c(kj2 + k^^) sin 



where k^ = k^ + k^ + k^^, c^ s [E(k) /2iTk^ ] , and ij; is a random number 
uniformly distributed on Che interval (0, 2tt). 


SCALING PROCEDURE 


The simulation variables are nondimensional (integer wave numbers with 
period 2 tt) and must be scaled to obtain dimensional values. If we wish to 
simulate an experimental flow, knowing at the initial time only its energy 
spectrum Ee(kg) and viscosity Vg» we must use a similar energy spectrum 
(i.e., differing only in energy and wave number scales) and specify a 
simulation viscosity such that the dimensionless problems are the same. Thus 
we define scale factors a, 3 relating the simulation energy E and wave 
number k to the experimental values by 

E(k)dk = aEg(k^)dkg , ke » 3k 

Here E and k are dimensionless, a has units L“^T^, and B has units 
L“^ so that time must scale as 
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and the* ^ 1 ^ ^ V 
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Thf scale factor u de))enda only on tlie scal ing, of the depend, nt. variables 
of the problem and. sine the computation allows an mvIlmUed ranp.e 
values cnani-ea of u produce only ctiani'ea of scale In the i eon Its. 1bl 
wl.uld also be the case for b M tlie computation allowed an unl.mited 
of values for the independent varial.les (wave runil)ers), but such would rcc.uin 

infinite spatial resolution. 

the computer simulation does not, of tourso, have Infinite tesolutlon, 
and the ranee of Its Independent variables is simply the value of Its l.lshesl 
“ondlmLsloLl (integer) Save number. Because the entire eaperlmenta range 
of length scales cannot be simulated, the choice of 6 determines mhto* 
physical wave numbers of the experiment are to be simulated. Clearly the 
error of the simulation depends on this choice, but at present no rationale 
is known for making it. The choice of Ferziger (ref. 1), 
the scaling that allows the greatest total energy to be included in the 
sLlatlo™ There Is also, however, tie Implicit 

period be "much greater" than the Integral scale of the turbulent flel . 

The number of modes, or degree of freedom of the motion, within a range 
of scales (e.g., K < k < 2K) characterized by K Is proportional to K . 

?L r»rereni^iale tanges of the motion ate certainly "ot represente egual y 
well in the statistical sense. There are very few modes in the ^ 

(smaller wave numbers), and if the simulation is to represent turbulence we 
must require that a small number of modes does not <:ontaln a large fractio 
of the energy. This is equivalent to the integral length constraint. 

If the truncation error is to be small, the viscosity must be high 
enough to damp the highest wave numbers of the calculation to the point where 
they, and presumably also those lost by truncation, have negligible effect. 

In other words, an accurate and complete solution requires that t p 

tion resolve all scales of motion; otherwise one must face the notoilous 
"closure problem." Since the simulation coJe presented here contains no 
closure approximations, It ie necessarily restricted to very low-turbulence 
Reynolds numbers. 

SAMPLE RESULTS 

Several runs of the simulation code have been mads for isotropic flow 
to develop the algorithm and to obtain numerical error estimates. 
cLrgy ere shown t„ figures 3-6. The algorithm used has negligible 

Lmerical dissipation and, when the spatial resolution In ..umated at auy 
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finite: wave: number In an Invlsclcl (v 0) fluid, 1 he- une'rgy cascade* eventually 
lt'>ads to an c'rjuip.irtil lone>d Mlaef* (fl|',. 1) with tlu* tbceeretlcal spectrum 
K(k) - ck^, 'liiJ:i ttni(kne:y Is at U 1 apuaMent wlie'ii iiioieculai* dissipation Is 
lncluele:d (flj’,. A) It the; range of computed scales does not Include almost 
all of the dlssli etion. More of the dissipation range can be Included 
by incre:aslng the range-* of computed scales (tig. 5) or by shifting the range 
of computed scales (fig. b) to Include more of tlic dissipation range. 
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THE ILL! AC PROGRAM 


Program Structure 

A fourth-order Runge-Kutta algorithia In used to integrate the system 
of equations (10)- (12). The strain inverses A, B, C, and the integrating 
factor F are considered known. The bulk of the computation is the evalua- 
tion of the right side of (10), whicti is^done in subroutines PHASEl, PHASE2, 
and PHASE3. The dependent variables X,Y,Z, are then advanced in STEP, and 
the continuity condition (11), is used by PRESSR to recover the physical 
velocities (12). These five subprograms are called sequentially by the 
control routine LOOP which is responsible for data management and step 
control. 

The functions of processes called by these routines are given by 
in-line comments in the listing. 


Data Structure and Flow 

The data base resides on disk and consists of two blocks. The first 
block of data holds the velocity field at the beginning of a Runge-Kutta 
step (three words/node) and a predicted velocity accumulator field 
(three words/node). This block of data is always accessed sequentially. 

The second block of data is working space (four words/node) in which the 
right side of (10) is evaluated, requiring both sequential and nonsequential 
page accesses from the disk. 

Each prediction within the Runge-Kutta process requires two complete 
passes through the data base, one bringing (x,y) planes into core (PHASEl, 
PHASES, STEP, and PRESSR) for operators in the y direction, and one bring- 
ing in (x,z) planes (PHASE2) for operators in the x and z directions. In 
the latter pass, only the working space data block is required, allowing 
the (x,z) planes to be handled by a triple buffered scheme. 


Listing of Program 

The program is coded for execution in 32-bit precision on the ILLIAC 
computer. The routines listed in this appendix, which are coded in the CFD 
language, cover the major algorithmic steps of the computation. Some of 
the lower level routines are coded in assembly language (ASK) for efficiency, 
and others had to be hand coded because of the restrictions placed on 32-blt 
operation by the CFD language. 
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Figure 2.- simulation program. 
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